** Clearing Stata memory
capture log close
clear all
set more off, perm
set seed 1234

///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////// Table O.47: Log (Annual Wages) between Seven and 12 Years After Admission Exam: Trimming /////
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

*****************************************************************************************************************************

** Opening Phase 2 norm_scores dataset 
use "Work Data/Gender_Phase2_long.dta",clear

*** Creating variables
encode subject, gen (sub)
tab subject, gen (d_sub)
label var sub "Subject"

** Subject dummies
rename d_sub1 Biology
rename d_sub2 Chemistry
rename d_sub3 Geography
rename d_sub4 History
rename d_sub5 Language
rename d_sub6 Mathematics
rename d_sub7 Physics
rename d_sub8 Portuguese
* Labels
label var Biology "Biology"
label var Chemistry "Chemistry"
label var Geography "Geography"
label var History "History"
label var Math "Mathematics"
label var Physics "Physics"
label var Portuguese "Portuguese"
label var Language "Foreign Language"

** Interaction: priority X female
gen fem_priority=female*priority
label var fem_priority "Female $\times$ Priority"

** Interaction: priority X subject
foreach v of varlist Biology-Portuguese {
gen fem_`v'=`v'*female
label var fem_`v' "Female $\times$ `v'"
gen prio_`v'=priority*`v'
label var prio_`v' "Priority $\times$ `v'"
gen fem_prio_`v'=fem_priority*`v'
label var fem_prio_`v' "Female $\times$ Priority $\times$ `v'"
}

global subject "Chemistry Geography History Mathematics Physics"
global subject_fem "fem_Chemistry fem_Geography fem_History fem_Mathematics fem_Physics"

** P1 scores: P1 normalized subject-specific scores
forvalues i=2(1)4 {
gen norm_p1score`i'=norm_p1score^`i'
sum norm_p1score`i'
}

*********************************************************************************
****************   Relative performances ****************************************
*********************************************************************************

******************************** ENEM ********************************

foreach v in norm_enem_w {
bys year female: egen `v'_ave_g=mean(`v')
gen `v'_g=`v'-`v'_ave_g
bys year female: sum `v'_g
}
drop norm_enem_w_ave_g

* Priority x relative performance in ENEM:
foreach v in norm_enem_w {
gen `v'_priority_g=`v'_g*priority
forvalues i=2(1)4 {
gen `v'_priority_g`i'=`v'_g^`i'*priority
sum `v'_priority_g`i'
}
}

global g_norm_enem_w_prio norm_enem_w_priority_g*
d $g_norm_enem_w_prio

* Interaction: subject X ENEM
foreach v of varlist Biology-Portuguese {
gen enem_`v'=`v'*norm_enem_w_g
label var enem_`v' "ENEM $\times$ `v'"
gen fem_enem_`v'=female*norm_enem_w_g*`v'
label var fem_enem_`v' "Female $\times$ ENEM $\times$ `v'"
forvalues i=2(1)4 {
gen enem_`v'_`i'=enem_`v'^`i'
gen fem_enem_`v'_`i'=fem_enem_`v'^`i'
}
sum enem_`v'* fem_enem_`v'*
}

global g_pol_enem_sub "enem_Chemistry* enem_Geography* enem_History* enem_Mathematics* enem_Physics*"
d $g_pol_enem_sub

******************************** Phase 1 scores ********************************

foreach v in norm_p1score {

tab year, sum(`v')
bys year subject female: egen gs_`v'_ave=mean(`v')
gen gs_`v'=`v'-gs_`v'_ave
bys year female subject: sum gs_`v'
drop gs_`v'_ave

forvalues i=2(1)4 {
gen gs_`v'`i'=gs_`v'^`i'
sum gs_`v'`i'
}

global gs_pol_`v' gs_`v' gs_`v'2 gs_`v'3 gs_`v'4
d $gs_pol_`v'

* Priority x Phase 1 scores:
gen gs_`v'_prio=gs_`v'*priority
forvalues i=2(1)4 {
gen gs_`v'_prio`i'=gs_`v'`i'*priority
sum gs_`v'_prio*
}
}

global gs_pol_norm_p1score_prio gs_norm_p1score_prio*
d $gs_pol_norm_p1score_prio

*********************************************************************************
**************** Main sample ****************************************************
*********************************************************************************

* 1) Only years before the affirmative action took place
drop if aa_year==1
tab year
drop if year==2000
tab year

* 2) Drop Portuguese and Foreign Language (in Phase 1 there is no Portuguese or Foreign Language exams - For Portuguese Phase 1 has an essay)
 tab subject, sum(norm_p1score)
 drop if subject=="lang" | subject=="port" 
 tab subject, sum(norm_p1score)
 drop Language Portuguese prio_Language prio_Portuguese fem_prio_Language fem_prio_Portuguese 
 
********************************************************************************************************
****************  Step 1: In the first step we obtain initial estimates and store them  ****************
********************************************************************************************************

** # 1 - Main specification, exclude interaction term
reghdfe norm_score priority $subject $subject_fem $g_pol_enem_sub $g_norm_enem_w_prio $gs_pol_norm_p1score $gs_pol_norm_p1score_prio, cluster(inscri2) absorb(inscri2) resid
** # 2 - Save residuals 
predict residuals, resid
** # 3 - Create measure of relative priority performance
ttest residuals, by(female)
ttest residuals if priority==1, by(female)
sum residuals if female==1 & priority==1
bys inscri2: egen res_prio_ave=mean(residuals) if priority==1
bys inscri2: egen res_prio_average=min(res_prio_ave)
mdesc res_prio_average
tab career_choice if  res_prio_average==.
bys inscri2: egen res_nonprio_ave=mean(residuals) if priority==0
bys inscri2: egen res_nonprio_average=min(res_nonprio_ave)
mdesc res_nonprio_average
gen res_diff=res_prio_average-res_nonprio_average
mdesc res_diff
ttest res_diff, by(female)

collapse (mean) res_diff female enem norm_enem_w gen_ques_st1 essay_st1 total_st1 career_choice year, by(inscri2)

count if res_diff==.

 ***************************************************
 ****************** Wages RAIS *********************
 ***************************************************

count
merge 1:1 inscri2 using "Work Data/RAIS_cleaned.dta"
drop if _merge==2
tab _merge  

** # 4 - Wage measures: 

forvalues i=7(1)12 {
gen yearafter`i'=year+`i'
tab yearafter`i', mi
gen annual_wage_after`i'=.
levelsof yearafter`i', local(levels) 
foreach l of local levels {
replace annual_wage_after`i'=mwagetot`l' if yearafter`i'==`l'
}
gen missing_wageafter`i' = missing(annual_wage_after`i')
tab missing_wageafter`i'
tab missing_wageafter`i' if _merge==1
}

forvalues i=7(1)12 {
count if missing_wageafter`i'==1 & december_wage_after`i'~=.
gen wagevalueafter`i'=1-missing_wageafter`i'
tab  wagevalueafter`i' missing_wageafter`i'
}

**** LFP between 7 and 12 years after admission exam

egen LFP712=rowmax(wagevalueafter7 wagevalueafter8 wagevalueafter9 wagevalueafter10 wagevalueafter11 wagevalueafter12) 
sum LFP712

**** Average and maximum wage

foreach x in annual_wage {
    
* Average
egen avg_712=rowmean(`x'_after7 `x'_after8 `x'_after9 `x'_after10 `x'_after11 `x'_after12) // 7-12 years after admission exam
sum avg_*

* Maximum
egen max_712=rowmax(`x'_after7 `x'_after8 `x'_after9 `x'_after10 `x'_after11 `x'_after12) // 7-12 years after admission exam
sum max_*

}

* Log variables
foreach v of varlist max* avg*  {
 gen l_`v'=log(`v')
 sum l_`v'
}

foreach v of varlist res_diff {
sum `v'
gen mean_`v'=r(mean)
gen sd_`v'=r(sd)
gen norm_`v'=(`v'-mean_`v')/sd_`v'
sum  norm_`v'
}

///////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////// BOUNDING /////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////

///////////////////////////////////////////////////////////////////////////////////////////
///////////////// 1) CHECK IMPACT ON SELECTION ///////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////

reghdfe LFP712 female norm_res_diff, vce(robust) absorb(year)
reghdfe LFP712 female norm_res_diff norm_enem_w, vce(robust) absorb(year)

** Control for program FE

reghdfe LFP712 female norm_res_diff, vce(robust) absorb(career_choice year)
reghdfe LFP712 female norm_res_diff norm_enem_w, vce(robust) absorb(career_choice year)

///////////////////////////////////////////////////////////////////////////////////////////
///////////////// 3) STORE BASELINE RESULTS //////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////

preserve

foreach y in LFP712 l_max_712 l_avg_712 { 	

reghdfe `y' female, vce(robust) absorb(year)
foreach x in female {
gen b`y'`x'_1=_b[`x']
gen se`y'`x'_1=_se[`x']
}

reghdfe `y' female norm_res_diff, vce(robust) absorb(year)
foreach x in female norm_res_diff {
gen b`y'`x'_2=_b[`x']
gen se`y'`x'_2=_se[`x']
}

reghdfe `y' female norm_enem_w, vce(robust) absorb(year)
foreach x in female norm_enem_w {
gen b`y'`x'_3=_b[`x']
gen se`y'`x'_3=_se[`x']
}

reghdfe `y' female norm_res_diff norm_enem_w, vce(robust) absorb(year)
foreach x in female norm_enem_w norm_res_diff {
gen b`y'`x'_4=_b[`x']
gen se`y'`x'_4=_se[`x']
}

reghdfe `y' female, vce(robust) absorb(career_choice year)
foreach x in female {
gen b`y'`x'_5=_b[`x']
gen se`y'`x'_5=_se[`x']
}

reghdfe `y' female norm_res_diff, vce(robust) absorb(career_choice year)
foreach x in female norm_res_diff {
gen b`y'`x'_6=_b[`x']
gen se`y'`x'_6=_se[`x']
}

reghdfe `y' female norm_enem_w, vce(robust) absorb(career_choice year)
foreach x in female norm_enem_w {
gen b`y'`x'_7=_b[`x']
gen se`y'`x'_7=_se[`x']
}

reghdfe `y' female norm_res_diff norm_enem_w, vce(robust) absorb(career_choice year)
foreach x in female norm_enem_w norm_res_diff {
gen b`y'`x'_8=_b[`x']
gen se`y'`x'_8=_se[`x']
}

}

keep b*_* se*_* 

duplicates drop

d

save "Output/baseline_results", replace

restore

///////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////// 2) BOUNDING ///////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////


*** Create variables for above and below mean relative priority performance
foreach x in norm_res_diff {
    egen median_`x'=median(`x')
	gen above_`x'=`x'>median_`x'
	sum median_`x' above_`x'
}

sum norm_res_diff if above_norm_res_diff==1

sum LFP712
gen mean_LFP_712=r(mean)

for any random norder: gen X=.


*************** Column (2), only control for the relative priority performance  ***************

* 1) Zero coefficient (trying to obtain coefficients lower than 0.001)

	*** Baseline effect
	foreach y in LFP712 l_max_712 l_avg_712 { 	
		reghdfe `y' female norm_res_diff, vce(robust) absorb(year)
		if "`y'"=="LFP712" {
			local beta = _b[norm_res_diff]
			local origbeta = _b[norm_res_diff]
		}	
	}
	
* Drop above (1) or below (0) median?
	if `origbeta'>0 local drop=1 // drop above median
	if `origbeta'<0 local drop=0 // drop below median
	di "`drop' and `origbeta'"

	*** Bounds
	set seed 12345

forvalues n=1(1)500 {
	
		preserve
		
		gen num=`n'
		local beta=`origbeta'
		
		* Create random reordoring of data so we drop at random
		replace random=runiform()
		sort above_norm_res_diff LFP712 random
		bysort above_norm_res_diff LFP712: replace norder=_n		
	
	local i=10
	count if above_norm_res_diff==1 
	
while (`i'<=3500 & (`beta'<-0.001 | `beta'>0.001)) & ((`origbeta'<0 & `beta'<0) | (`origbeta'>0 & `beta'>0))	{
replace l_max_712=. if LFP712!=0 & above_norm_res_diff==`drop' & norder<`i'
replace l_avg_712=. if LFP712!=0 & above_norm_res_diff==`drop' & norder<`i'
replace LFP712=0 if LFP712!=0  & above_norm_res_diff==`drop' & norder<`i'
quietly	reghdfe LFP712 female norm_res_diff, vce(robust) absorb(year)
			local beta = _b[norm_res_diff]
			di "Beta interaction on LFP for dropping `i' is `beta'"

	local i=`i'+10	
}

* Save new coefficients
foreach y in  LFP712 l_max_712 l_avg_712 { 	
quietly	reghdfe `y' female , vce(robust) absorb(year)
foreach x in female {
gen b`y'`x'_1=_b[`x']
}
quietly	reghdfe `y' female norm_res_diff, vce(robust) absorb(year)
	foreach x in female norm_res_diff {
gen b`y'`x'_2=_b[`x']
}
			if "`y'"=="l_avg_712" gen wage_obs=e(N)
		}
		collapse (mean) b*_* wage_obs, by(num)
		sum		
		if num==1 save Output/rand_bounds_col2, replace
		if num>1 {
			append using Output/rand_bounds_col2
			sleep 50
			save Output/rand_bounds_col2, replace
		}
		
		restore

	}	

*************** Column (4),  control for the relative priority performance and ENEM scores ***************	

* 1) Zero coefficient (trying to obtain coefficients lower than 0.001)

	*** Baseline effect
	foreach y in LFP712 l_max_712 l_avg_712 { 	
		reghdfe `y' female norm_res_diff norm_enem_w, vce(robust) absorb(year)
		if "`y'"=="LFP712" {
			local beta = _b[norm_res_diff]
			local origbeta = _b[norm_res_diff]
		}	
	}
	
* Drop above (1) or below (0) median?
	if `origbeta'>0 local drop=1 // drop above median
	if `origbeta'<0 local drop=0 // drop below median
	di "`drop' and `origbeta'"

	*** Bounds
	set seed 12345

forvalues n=1(1)500 {
	
		preserve
		
		gen num=`n'
		local beta=`origbeta'
		
		* Create random reording of data so we drop at random
		replace random=runiform()
		sort above_norm_res_diff LFP712 random
		bysort above_norm_res_diff LFP712: replace norder=_n		
	
	local i=10
	count if above_norm_res_diff==1 
	
while (`i'<=3500 & (`beta'<-0.001 | `beta'>0.001)) & ((`origbeta'<0 & `beta'<0) | (`origbeta'>0 & `beta'>0))	{
replace l_max_712=. if LFP712!=0 & above_norm_res_diff==`drop' & norder<`i'
replace l_avg_712=. if LFP712!=0 & above_norm_res_diff==`drop' & norder<`i'
replace LFP712=0 if LFP712!=0  & above_norm_res_diff==`drop' & norder<`i'
quietly	reghdfe LFP712 female norm_res_diff norm_enem_w, vce(robust) absorb(year)
			local beta = _b[norm_res_diff]
			di "Beta interaction on LFP for dropping `i' is `beta'"

	local i=`i'+10	
}
* Save new coefficients
foreach y in  LFP712 l_max_712 l_avg_712 { 	

quietly	reghdfe `y' female norm_enem_w, vce(robust) absorb(year)
foreach x in female norm_enem_w {
gen b`y'`x'_3=_b[`x']
}
			
quietly	reghdfe `y' female norm_res_diff norm_enem_w, vce(robust) absorb(year)
foreach x in female norm_enem_w norm_res_diff {
gen b`y'`x'_4=_b[`x']
}	
			if "`y'"=="l_avg_712" gen wage_obs=e(N)
		}
		collapse (mean) b*_* wage_obs, by(num)
		sum		
		if num==1 save Output/rand_bounds_col4, replace
		if num>1 {
			append using Output/rand_bounds_col4
			sleep 50
			save Output/rand_bounds_col4, replace
		}
		
		restore

	}	
	

*************** Column (6),  control for the relative priority performance and major FE ***************	

* 1) Zero coefficient (trying to obtain coefficients lower than 0.001)

	*** Baseline effect
	foreach y in LFP712 l_max_712 l_avg_712 { 	
		reghdfe `y' female norm_res_diff, vce(robust) absorb(career_choice year)
		if "`y'"=="LFP712" {
			local beta = _b[norm_res_diff]
			local origbeta = _b[norm_res_diff]
		}	
	}
	
* Drop above (1) or below (0) median?
	if `origbeta'>0 local drop=1 // drop above median
	if `origbeta'<0 local drop=0 // drop below median
	di "`drop' and `origbeta'"

	*** Bounds
	set seed 12345

forvalues n=1(1)500 {
	
		preserve
		
		gen num=`n'
		local beta=`origbeta'
		
		* Create random reording of data so we drop at random
		replace random=runiform()
		sort above_norm_res_diff LFP712 random
		bysort above_norm_res_diff LFP712: replace norder=_n		
	
	local i=10
	count if above_norm_res_diff==1 
	
while (`i'<=3500 & (`beta'<-0.001 | `beta'>0.001)) & ((`origbeta'<0 & `beta'<0) | (`origbeta'>0 & `beta'>0))	{
replace l_max_712=. if LFP712!=0 & above_norm_res_diff==`drop' & norder<`i'
replace l_avg_712=. if LFP712!=0 & above_norm_res_diff==`drop' & norder<`i'
replace LFP712=0 if LFP712!=0  & above_norm_res_diff==`drop' & norder<`i'
quietly	reghdfe LFP712 female norm_res_diff, vce(robust) absorb(career_choice year)
			local beta = _b[norm_res_diff]
			di "Beta interaction on LFP for dropping `i' is `beta'"

	local i=`i'+10	
}
* Save new coefficients
foreach y in  LFP712 l_max_712 l_avg_712 { 	
	
quietly	reghdfe `y' female , vce(robust) absorb(career_choice year)
foreach x in female {
gen b`y'`x'_5=_b[`x']
}
		
quietly	reghdfe `y' female norm_res_diff, vce(robust) absorb(career_choice year)
foreach x in female norm_res_diff {
gen b`y'`x'_6=_b[`x']
}
if "`y'"=="l_avg_712" gen wage_obs=e(N)

}
		collapse (mean) b*_* wage_obs, by(num)
		sum		
		if num==1 save Output/rand_bounds_col6, replace
		if num>1 {
			append using Output/rand_bounds_col6
			sleep 50
			save Output/rand_bounds_col6, replace
		}
		
		restore

	}	


*************** Column (8),  control for the relative priority performance, ENEM and major FE ***************	

* 1) Zero coefficient (trying to obtain coefficients lower than 0.001)

	*** Baseline effect
	foreach y in LFP712 l_max_712 l_avg_712 { 	
		reghdfe `y' female norm_res_diff norm_enem_w, vce(robust) absorb(career_choice year)
		if "`y'"=="LFP712" {
			local beta = _b[norm_res_diff]
			local origbeta = _b[norm_res_diff]
		}	
	}
	
* Drop above (1) or below (0) median?
	if `origbeta'>0 local drop=1 // drop above median
	if `origbeta'<0 local drop=0 // drop below median
	di "`drop' and `origbeta'"

	*** Bounds
	set seed 12345

forvalues n=1(1)500 {
	
		preserve
		
		gen num=`n'
		local beta=`origbeta'
		
		* Create random reording of data so we drop at random
		replace random=runiform()
		sort above_norm_res_diff LFP712 random
		bysort above_norm_res_diff LFP712: replace norder=_n		
	
	local i=10
	count if above_norm_res_diff==1 
	
while (`i'<=3500 & (`beta'<-0.001 | `beta'>0.001)) & ((`origbeta'<0 & `beta'<0) | (`origbeta'>0 & `beta'>0))	{
replace l_max_712=. if LFP712!=0 & above_norm_res_diff==`drop' & norder<`i'
replace l_avg_712=. if LFP712!=0 & above_norm_res_diff==`drop' & norder<`i'
replace LFP712=0 if LFP712!=0  & above_norm_res_diff==`drop' & norder<`i'
quietly	reghdfe LFP712 female norm_res_diff norm_enem_w, vce(robust) absorb(career_choice year)
			local beta = _b[norm_res_diff]
			di "Beta interaction on LFP for dropping `i' is `beta'"

	local i=`i'+10	
}
* Save new coefficients
foreach y in  LFP712 l_max_712 l_avg_712 { 	
	
quietly	reghdfe `y' female  norm_enem_w, vce(robust) absorb(career_choice year)
foreach x in female norm_enem_w {
gen b`y'`x'_7=_b[`x']
}
			
quietly	reghdfe `y' female norm_res_diff norm_enem_w, vce(robust) absorb(career_choice year)
foreach x in female norm_enem_w norm_res_diff {
gen b`y'`x'_8=_b[`x']
}
			
			if "`y'"=="l_avg_712" gen wage_obs=e(N)
		}
		collapse (mean) b*_* wage_obs, by(num)
		sum		
		if num==1 save Output/rand_bounds_col8, replace
		if num>1 {
			append using Output/rand_bounds_col8
			sleep 50
			save Output/rand_bounds_col8, replace
		}
		
		restore

	}	
	
** Append data files from different columns
use Output/rand_bounds_col2, clear
reshape long bLFP712female_ bLFP712norm_res_diff_ bl_max_712female_  bl_max_712norm_res_diff_ bl_avg_712female_ bl_avg_712norm_res_diff_, i(num) j(col)
gen bound="col12"
save Output/rand_bounds_col2_long

use Output/rand_bounds_col4, clear
reshape long bLFP712female_ bLFP712norm_res_diff_ bLFP712norm_enem_w_ bl_max_712female_  bl_max_712norm_res_diff_  bl_max_712norm_enem_w_ bl_avg_712female_ bl_avg_712norm_res_diff_   bl_avg_712norm_enem_w_, i(num) j(col)
gen bound="col34"
save Output/rand_bounds_col4_long

use Output/rand_bounds_col6, clear
reshape long bLFP712female_ bLFP712norm_res_diff_ bl_max_712female_  bl_max_712norm_res_diff_ bl_avg_712female_ bl_avg_712norm_res_diff_, i(num) j(col)
gen bound="col56"
save Output/rand_bounds_col6_long

use Output/rand_bounds_col8, clear
reshape long bLFP712female_ bLFP712norm_res_diff_ bLFP712norm_enem_w_ bl_max_712female_  bl_max_712norm_res_diff_  bl_max_712norm_enem_w_ bl_avg_712female_ bl_avg_712norm_res_diff_   bl_avg_712norm_enem_w_, i(num) j(col)
gen bound="col78"
save Output/rand_bounds_col8_long

use Output/rand_bounds_col2_long
forvalues i=4(2)8{
append using Output/rand_bounds_col`i'_long
tab col
}
tab col
save Output/rand_bounds, replace


************ Export table ************

use Output/rand_bounds, clear

* Formal labor market participation
 
 * 95% confidence intervals
 
label var bl_avg_712female_ "Female"
label var bl_max_712female_ "Female"
label var bl_avg_712norm_res_diff_ "Relative priority performance"
label var bl_avg_712norm_enem_w_ "Norm. ENEM scores"
label var bl_max_712norm_res_diff_ "Relative priority performance"
label var bl_max_712norm_enem_w_ "Norm. ENEM scores"

* Avg
 
estpost sum bl_avg_712female_  if col==1
eststo col1
centile bl_avg_712female_  if col==1, centile(2.5 97.5)
matrix lb = r(c_1)
matrix ub = r(c_2)
mat colnames lb = "bl_avg_712female_"
mat colnames ub = "bl_avg_712female_"
estadd matrix lb = lb
estadd matrix ub = ub

estpost sum bl_avg_712female_ bl_avg_712norm_res_diff_  if col==2
eststo col2
centile bl_avg_712female_ if col==2, centile(2.5 97.5)
 matrix lb1 = r(c_1)
 matrix ub1 = r(c_2)
centile bl_avg_712norm_res_diff_ if col==2, centile(2.5 97.5)
 matrix lb2 = r(c_1)
 matrix ub2 = r(c_2)
 matrix lb = lb1 , lb2
 matrix ub = ub1 , ub2
mat colnames lb = "bl_avg_712female_" "bl_avg_712norm_res_diff_"
mat colnames ub = "bl_avg_712female_" "bl_avg_712norm_res_diff_"
estadd matrix lb = lb
estadd matrix ub = ub
 
estpost sum bl_avg_712female_ bl_avg_712norm_enem_w_ if col==3
eststo col3
centile bl_avg_712female_ if col==3, centile(2.5 97.5)
 matrix lb1 = r(c_1)
 matrix ub1 = r(c_2)
centile bl_avg_712norm_enem_w_ if col==3, centile(2.5 97.5)
 matrix lb2 = r(c_1)
 matrix ub2 = r(c_2)
 matrix lb = lb1 , lb2
 matrix ub = ub1 , ub2
 mat colnames lb = "bl_avg_712female_" "bl_avg_712norm_enem_w_"
 mat colnames ub = "bl_avg_712female_" "bl_avg_712norm_enem_w_"
estadd matrix lb = lb
estadd matrix ub = ub

estpost sum bl_avg_712female_ bl_avg_712norm_res_diff_ bl_avg_712norm_enem_w_ if col==4
eststo col4
centile bl_avg_712female_  if col==4, centile(2.5 97.5)
 matrix lb1 = r(c_1)
 matrix ub1 = r(c_2)
centile bl_avg_712norm_res_diff_ if col==4, centile(2.5 97.5)
 matrix lb2 = r(c_1)
 matrix ub2 = r(c_2)
 centile bl_avg_712norm_enem_w_ if col==4, centile(2.5 97.5)
 matrix lb3= r(c_1)
 matrix ub3 = r(c_2)
 matrix lb = lb1 , lb2 , lb3
 matrix ub = ub1 , ub2 , ub3
 mat colnames lb = "bl_avg_712female_" "bl_avg_712norm_res_diff_" "bl_avg_712norm_enem_w_"
 mat colnames ub = "bl_avg_712female_" "bl_avg_712norm_res_diff_" "bl_avg_712norm_enem_w_"
estadd matrix lb = lb
estadd matrix ub = ub

estpost sum bl_avg_712female_ if col==5
eststo col5
centile bl_avg_712female_  if col==5, centile(2.5 97.5)
matrix lb = r(c_1)
matrix ub = r(c_2)
mat colnames lb = "bl_avg_712female_"
mat colnames ub = "bl_avg_712female_"
estadd matrix lb = lb
estadd matrix ub = ub

estpost sum bl_avg_712female_ bl_avg_712norm_res_diff_  if col==6
eststo col6
centile bl_avg_712female_  if col==6, centile(2.5 97.5)
 matrix lb1 = r(c_1)
 matrix ub1 = r(c_2)
centile bl_avg_712norm_res_diff_ if col==6, centile(2.5 97.5)
 matrix lb2 = r(c_1)
 matrix ub2 = r(c_2)
 matrix lb = lb1 , lb2
 matrix ub = ub1 , ub2
mat colnames lb = "bl_avg_712female_" "bl_avg_712norm_res_diff_"
mat colnames ub = "bl_avg_712female_" "bl_avg_712norm_res_diff_"
estadd matrix lb = lb
estadd matrix ub = ub

estpost sum bl_avg_712female_ bl_avg_712norm_enem_w_ if col==7
eststo col7
centile bl_avg_712female_ if col==7, centile(2.5 97.5)
 matrix lb1 = r(c_1)
 matrix ub1 = r(c_2)
centile bl_avg_712norm_enem_w_ if col==7, centile(2.5 97.5)
 matrix lb2 = r(c_1)
 matrix ub2 = r(c_2)
 matrix lb = lb1 , lb2
 matrix ub = ub1 , ub2
 mat colnames lb = "bl_avg_712female_" "bl_avg_712norm_enem_w_"
 mat colnames ub = "bl_avg_712female_" "bl_avg_712norm_enem_w_"
estadd matrix lb = lb
estadd matrix ub = ub

estpost sum bl_avg_712female_ bl_avg_712norm_res_diff_ bl_avg_712norm_enem_w_ if col==8
eststo col8
centile bl_avg_712female_ if col==8, centile(2.5 97.5)
 matrix lb1 = r(c_1)
 matrix ub1 = r(c_2)
centile bl_avg_712norm_res_diff_ if col==8, centile(2.5 97.5)
 matrix lb2 = r(c_1)
 matrix ub2 = r(c_2)
 centile bl_avg_712norm_enem_w_ if col==8, centile(2.5 97.5)
 matrix lb3= r(c_1)
 matrix ub3 = r(c_2)
 matrix lb = lb1 , lb2 , lb3
 matrix ub = ub1 , ub2 , ub3
 mat colnames lb = "bl_avg_712female_" "bl_avg_712norm_res_diff_" "bl_avg_712norm_enem_w_"
 mat colnames ub = "bl_avg_712female_" "bl_avg_712norm_res_diff_" "bl_avg_712norm_enem_w_"
estadd matrix lb = lb
estadd matrix ub = ub

esttab col* using "Output\Bound_Wages.tex", cells(mean(fmt(3)) lb(par("[") fmt(3))&ub(fmt(3))&count(par("]") fmt(0))) incelldelimiter(",") substitute (",]500" "]" ",," "") label collabels(none) replace noobs booktabs nomtitle nogaps f refcat(bl_avg_712female_ " \\ \multicolumn{9}{l}{\textit{Panel A: Average annual wages (7-12 years after admission exam)}} \\", nolabel)

* Maximum

estpost sum bl_max_712female_  if col==1
eststo col1
centile bl_max_712female_  if col==1, centile(2.5 97.5)
matrix lb = r(c_1)
matrix ub = r(c_2)
mat colnames lb = "bl_max_712female_"
mat colnames ub = "bl_max_712female_"
estadd matrix lb = lb
estadd matrix ub = ub

estpost sum bl_max_712female_ bl_max_712norm_res_diff_  if col==2
eststo col2
centile bl_max_712female_ if col==2, centile(2.5 97.5)
 matrix lb1 = r(c_1)
 matrix ub1 = r(c_2)
centile bl_max_712norm_res_diff_ if col==2, centile(2.5 97.5)
 matrix lb2 = r(c_1)
 matrix ub2 = r(c_2)
 matrix lb = lb1 , lb2
 matrix ub = ub1 , ub2
mat colnames lb = "bl_max_712female_" "bl_max_712norm_res_diff_"
mat colnames ub = "bl_max_712female_" "bl_max_712norm_res_diff_"
estadd matrix lb = lb
estadd matrix ub = ub
 
estpost sum bl_max_712female_ bl_max_712norm_enem_w_ if col==3
eststo col3
centile bl_max_712female_ if col==3, centile(2.5 97.5)
 matrix lb1 = r(c_1)
 matrix ub1 = r(c_2)
centile bl_max_712norm_enem_w_ if col==3, centile(2.5 97.5)
 matrix lb2 = r(c_1)
 matrix ub2 = r(c_2)
 matrix lb = lb1 , lb2
 matrix ub = ub1 , ub2
 mat colnames lb = "bl_max_712female_" "bl_max_712norm_enem_w_"
 mat colnames ub = "bl_max_712female_" "bl_max_712norm_enem_w_"
estadd matrix lb = lb
estadd matrix ub = ub

estpost sum bl_max_712female_ bl_max_712norm_res_diff_ bl_max_712norm_enem_w_ if col==4
eststo col4
centile bl_max_712female_  if col==4, centile(2.5 97.5)
 matrix lb1 = r(c_1)
 matrix ub1 = r(c_2)
centile bl_max_712norm_res_diff_ if col==4, centile(2.5 97.5)
 matrix lb2 = r(c_1)
 matrix ub2 = r(c_2)
 centile bl_max_712norm_enem_w_ if col==4, centile(2.5 97.5)
 matrix lb3= r(c_1)
 matrix ub3 = r(c_2)
 matrix lb = lb1 , lb2 , lb3
 matrix ub = ub1 , ub2 , ub3
 mat colnames lb = "bl_max_712female_" "bl_max_712norm_res_diff_" "bl_max_712norm_enem_w_"
 mat colnames ub = "bl_max_712female_" "bl_max_712norm_res_diff_" "bl_max_712norm_enem_w_"
estadd matrix lb = lb
estadd matrix ub = ub

estpost sum bl_max_712female_ if col==5
eststo col5
centile bl_max_712female_  if col==5, centile(2.5 97.5)
matrix lb = r(c_1)
matrix ub = r(c_2)
mat colnames lb = "bl_max_712female_"
mat colnames ub = "bl_max_712female_"
estadd matrix lb = lb
estadd matrix ub = ub

estpost sum bl_max_712female_ bl_max_712norm_res_diff_  if col==6
eststo col6
centile bl_max_712female_  if col==6, centile(2.5 97.5)
 matrix lb1 = r(c_1)
 matrix ub1 = r(c_2)
centile bl_max_712norm_res_diff_ if col==6, centile(2.5 97.5)
 matrix lb2 = r(c_1)
 matrix ub2 = r(c_2)
 matrix lb = lb1 , lb2
 matrix ub = ub1 , ub2
mat colnames lb = "bl_max_712female_" "bl_max_712norm_res_diff_"
mat colnames ub = "bl_max_712female_" "bl_max_712norm_res_diff_"
estadd matrix lb = lb
estadd matrix ub = ub

estpost sum bl_max_712female_ bl_max_712norm_enem_w_ if col==7
eststo col7
centile bl_max_712female_ if col==7, centile(2.5 97.5)
 matrix lb1 = r(c_1)
 matrix ub1 = r(c_2)
centile bl_max_712norm_enem_w_ if col==7, centile(2.5 97.5)
 matrix lb2 = r(c_1)
 matrix ub2 = r(c_2)
 matrix lb = lb1 , lb2
 matrix ub = ub1 , ub2
 mat colnames lb = "bl_max_712female_" "bl_max_712norm_enem_w_"
 mat colnames ub = "bl_max_712female_" "bl_max_712norm_enem_w_"
estadd matrix lb = lb
estadd matrix ub = ub

estpost sum bl_max_712female_ bl_max_712norm_res_diff_ bl_max_712norm_enem_w_ if col==8
eststo col8
centile bl_max_712female_ if col==8, centile(2.5 97.5)
 matrix lb1 = r(c_1)
 matrix ub1 = r(c_2)
centile bl_max_712norm_res_diff_ if col==8, centile(2.5 97.5)
 matrix lb2 = r(c_1)
 matrix ub2 = r(c_2)
 centile bl_max_712norm_enem_w_ if col==8, centile(2.5 97.5)
 matrix lb3= r(c_1)
 matrix ub3 = r(c_2)
 matrix lb = lb1 , lb2 , lb3
 matrix ub = ub1 , ub2 , ub3
 mat colnames lb = "bl_max_712female_" "bl_max_712norm_res_diff_" "bl_max_712norm_enem_w_"
 mat colnames ub = "bl_max_712female_" "bl_max_712norm_res_diff_" "bl_max_712norm_enem_w_"
estadd matrix lb = lb
estadd matrix ub = ub

estadd local program "Yes":col5 col6 col7 col8
estadd local program "No":col1 col2 col3 col4
estadd local year "Yes": col*

 esttab col* using "Output\Bound_Wages.tex", cells(mean(fmt(3)) lb(par("[") fmt(3))&ub(fmt(3))&count(par("]") fmt(0))) incelldelimiter(",") substitute (",]500" "]" ",," "") label collabels(none) append stats(N year program, fmt( %9.0fc %3s %3s) labels("Number of random draws" "Exam year FE" "Major FE")) booktabs nomtitle nogaps f refcat(bl_max_712female_ " \\ \multicolumn{9}{l}{\textit{Panel B: Maximum annual wages (7-12 years after admission exam)}} \\", nolabel) nonumber
 